PHC 6937 (Data Visualization) Final Project

Introduction

Libaries

library(dplyr)       # tidy data
library(Seurat)      # scRNA-seq infrastructure
library(ggpubr)      # plot utilities
library(celldex)     # RNA-seq reference datasets
library(ggplot2)     # plots
library(SingleR)     # cell annotations
library(decoderr)    # pathway analysis through NMF
library(SCISSORS)    # recluster scRNA-seq data
library(CONICSmat)   # estimate CNVs
library(patchwork)   # arrange plots
library(jackknife)   # my package
library(paletteer)   # all the palettes
library(kableExtra)  # table styles

Data

There are several samples here, so we’ll need to create a dictionary of the samplenames and the sample types.

dir_names <- c("GSM4679532", "GSM4679533", "GSM4679534", "GSM4679535", "GSM4679536", 
               "GSM4679537", "GSM4679538", "GSM4679539", "GSM4679540", "GSM4679541", 
               "GSM4679542", "GSM4679543", "GSM4679544", "GSM4679545", "GSM4679546", "GSM4679547")
sample_names <- c("P01", "P02", "P03", "P04", "P05", "P06", "P07", "P08", "P09", "P10", 
                  "MET01", "MET02", "MET03", "MET04", "MET05", "MET06")
sample_types <- c(rep("Primary", 10), rep("Metastatic", 6))
sample_dict <- data.frame(Directory = dir_names, 
                          Sample = sample_names, 
                          Type = sample_types)
sample_dict %>% 
  kbl(booktabs = TRUE, caption = "Sample Dictionary") %>% 
  kable_minimal("hover", html_font = "Avenir", full_width = FALSE)

Sample Dictionary

Directory	Sample	Type
GSM4679532	P01	Primary
GSM4679533	P02	Primary
GSM4679534	P03	Primary
GSM4679535	P04	Primary
GSM4679536	P05	Primary
GSM4679537	P06	Primary
GSM4679538	P07	Primary
GSM4679539	P08	Primary
GSM4679540	P09	Primary
GSM4679541	P10	Primary
GSM4679542	MET01	Metastatic
GSM4679543	MET02	Metastatic
GSM4679544	MET03	Metastatic
GSM4679545	MET04	Metastatic
GSM4679546	MET05	Metastatic
GSM4679547	MET06	Metastatic

Dataset Integration

First, we loop across the samples, create Seurat objects, and normalize expression and identify variable genes. The first 10 samples are primary tumors, and the next 6 are metastatic.

obj_list <- list()
for (i in seq(sample_dict$Directory)) {
  file <- paste0("~/Desktop/Data/Lin et al/GSE154778_RAW/", sample_dict$Directory[i])
  counts <- Read10X(file)
  seurat_obj <- CreateSeuratObject(counts, 
                                   project = sample_dict$Sample[i], 
                                   min.cells = 3, 
                                   min.features = 1000)
  seurat_obj@meta.data$sample <- sample_dict$Sample[i]
  seurat_obj[["percent_MT"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-")
  seurat_obj <- SCTransform(seurat_obj, 
                            variable.features.n = 4000, 
                            vars.to.regress = "percent_MT", 
                            seed.use = 629, 
                            verbose = FALSE)
  obj_list[[i]] <- seurat_obj
}
primary_obj <- obj_list[1:10]
meta_obj <- obj_list[11:16]

Next we integrate all the samples into one Seurat object. We’ll start with the primary samples.

integration_genes <- SelectIntegrationFeatures(object.list = primary_obj, 
                                               nfeatures = 4000, 
                                               verbose = FALSE)
primary_obj <- PrepSCTIntegration(object.list = primary_obj, 
                                  anchor.features = integration_genes, 
                                  verbose = FALSE)
integration_anchors <- FindIntegrationAnchors(object.list = primary_obj, 
                                              normalization.method = "SCT", 
                                              anchor.features = integration_genes, 
                                              verbose = FALSE)
primary <- IntegrateData(anchorset = integration_anchors, 
                         normalization.method = "SCT", 
                         k.weight = 50,
                         verbose = FALSE)
primary <- CellCycleScoring(primary, 
                            s.features = cc.genes.updated.2019$s.genes, 
                            g2m.features = cc.genes.updated.2019$g2m.genes)
DefaultAssay(primary) <- "integrated"

We repeat the process for the metastatic samples.

integration_genes <- SelectIntegrationFeatures(object.list = meta_obj, 
                                               nfeatures = 4000, 
                                               verbose = FALSE)
meta_obj <- PrepSCTIntegration(object.list = meta_obj, 
                               anchor.features = integration_genes, 
                               verbose = FALSE)
integration_anchors <- FindIntegrationAnchors(object.list = meta_obj, 
                                              normalization.method = "SCT", 
                                              anchor.features = integration_genes, 
                                              verbose = FALSE)
meta <- IntegrateData(anchorset = integration_anchors, 
                      normalization.method = "SCT", 
                      k.weight = 50,
                      verbose = FALSE)
meta <- CellCycleScoring(meta, 
                         s.features = cc.genes.updated.2019$s.genes, 
                         g2m.features = cc.genes.updated.2019$g2m.genes)
DefaultAssay(meta) <- "integrated"

Preprocessing

Primary Tumors

We run PCA / UMAP, create the SNN graph representation of our cells, and cluster the cells using Louvain modularity optimization.

primary <- RunPCA(primary, 
                  npcs = 30, 
                  verbose = FALSE, 
                  seed.use = 629)
primary <- RunUMAP(primary, 
                   reduction = "pca", 
                   dims = 1:15, 
                   metric = "cosine",
                   min.dist = 0.1, 
                   n.components = 2, 
                   seed.use = 629, 
                   verbose = FALSE)
primary <- FindNeighbors(primary, 
                         reduction = "pca", 
                         dims = 1:15, 
                         k.param = 30, 
                         annoy.metric = "cosine", 
                         verbose = FALSE)
primary <- FindClusters(primary, 
                        resolution = 0.3, 
                        random.seed = 629,
                        verbose = FALSE)

Metastatic Tumors

Now we repeat the preprocessing on the metastatic tumors.

meta <- RunPCA(meta, 
               npcs = 30, 
               verbose = FALSE, 
               seed.use = 629)
meta <- RunUMAP(meta, 
                reduction = "pca", 
                dims = 1:30, 
                min.dist = 0.05,
                metric = "cosine", 
                n.components = 2, 
                seed.use = 629, 
                verbose = FALSE)
meta <- FindNeighbors(meta, 
                      reduction = "pca", 
                      dims = 1:30, 
                      annoy.metric = "cosine", 
                      verbose = FALSE)
meta <- FindClusters(meta, 
                     resolution = 0.15, 
                     verbose = FALSE, 
                     random.seed = 629)

Visualization

First we visualize the semi-supervised clustering we’ve generated using Louvain modularity optimization.

p0 <- DimPlot(primary) + 
      labs(subtitle = "Primary Tumor", x = "UMAP 1", y = "UMAP 2") + 
      scale_color_paletteer_d("ggthemes::Classic_20") + 
      theme_yehlab() + 
      theme(axis.title.x = element_text(hjust = 0, size = 10), 
            axis.title.y = element_text(hjust = 0, size = 10), 
            plot.subtitle = element_text(hjust = 0.5, size = 12)) + 
      guides(color = guide_legend(nrow = 2, byrow = TRUE, override.aes = list(size = 3)))
p1 <- DimPlot(meta, cols = paletteer_d("ggthemes::Classic_20")[8:12]) + 
      labs(subtitle = "Metastatic Tumor") + 
      theme_yehlab() + 
      theme(axis.title = element_blank(), 
            plot.subtitle = element_text(hjust = 0.5, size = 12)) + 
      guides(color = guide_legend(nrow = 1, byrow = TRUE, override.aes = list(size = 3)))
(p0 | p1) + plot_annotation(title = "Louvain Clustering", 
                            theme = theme(plot.title = element_text(hjust = 0.5)))

Next we check out how the samples are arranged in UMAP space. We see that there is next to no grouping by sample, which tells us that the data integration process worked correctly and the sample-related batch effect has been regressed out.

p2 <- DimPlot(primary, group.by = "sample") + 
      labs(x = "UMAP 1", y = "UMAP 2", subtitle = "Primary Tumor") + 
      scale_color_paletteer_d("miscpalettes::pastel") + 
      theme_yehlab() + 
      theme(axis.title = element_text(hjust = 0, size = 10), 
            plot.subtitle = element_text(hjust = 0.5), 
            legend.text = element_text(size = 10), 
            plot.title = element_blank()) + 
      guides(color = guide_legend(nrow = 2, byrow = TRUE, override.aes = list(size = 3)))
p3 <- DimPlot(meta, group.by = "sample") + 
      labs(subtitle = "Metastatic Tumor") + 
      scale_color_paletteer_d("miscpalettes::pastel") + 
      theme_yehlab() + 
      theme(axis.title = element_blank(), 
            plot.subtitle = element_text(hjust = 0.5), 
            legend.text = element_text(size = 10), 
            plot.title = element_blank()) + 
      guides(color = guide_legend(nrow = 2, byrow = TRUE, override.aes = list(size = 3)))
(p2 | p3) + plot_annotation(title = "Sample IDs", 
                            theme = theme(plot.caption = element_text(face = "italic"), 
                                          plot.title = element_text(hjust= 0.5)))

Analysis

Differential Expression

Primary Tumor

We use a Wilcox test for differential expression, and filter the results after the Bonferroni p-value correction at the \(q < 0.01\) level.

primary_markers <- FindAllMarkers(primary, 
                                  logfc.threshold = 2, 
                                  test.use = "wilcox", 
                                  only.pos = TRUE, 
                                  random.seed = 629, 
                                  verbose = FALSE) %>% 
                   filter(p_val_adj < 0.01)
primary_markers %>% 
  group_by(cluster) %>% 
  top_n(n = 5, wt = avg_log2FC) -> primary_markers_top5
p4 <- DotPlot(primary, 
              cols = c("grey70", "firebrick"), 
              features = primary_markers_top5$gene) + 
      labs(x = "Gene", 
           y = "Cluster", 
           color = "Mean Expression", 
           size = "Percent Expressed", 
           title = "Primary Tumor Marker Genes") + 
      theme(legend.position = "bottom", 
            axis.line = element_blank(), 
            legend.justification = "center", 
            plot.title = element_text(hjust = 0.5), 
            axis.title.y = element_text(angle = 360, vjust = 0.5), 
            panel.border = element_rect(fill = NA, size = 1, color = "black"), 
            axis.text.x = element_text(angle = 45, size = 8, vjust = 1, hjust = 1)) + 
      guides(color = guide_colorbar(title.position = "top", title.hjust = 0.5, barwidth = unit(5, "cm")), 
             size = guide_legend(title.position = "top", title.hjust = 0.5))

Clusters 0 and 1 are partially determined by LYZ & CXCL5, respectively, both of which are markers for Epithelial cells. Cluster 3 has CD68, which is a myeloid marker.

p4

Metastatic Tumor

We perform the same DE test on the metastatic tumor cells.

meta_markers <- FindAllMarkers(meta, 
                               logfc.threshold = 2, 
                               test.use = "wilcox", 
                               only.pos = TRUE, 
                               random.seed = 629, 
                               verbose = FALSE) %>% 
                filter(p_val_adj < 0.01)
meta_markers %>% 
  group_by(cluster) %>% 
  top_n(n = 5, wt = avg_log2FC) -> meta_markers_top5
p5 <- DotPlot(meta, 
              cols = c("grey70", "firebrick"), 
              features = meta_markers_top5$gene) + 
      labs(x = "Gene", 
           y = "Cluster", 
           color = "Mean Expression", 
           size = "Percent Expressed", 
           title = "Metastatic Tumor Marker Genes") + 
      theme(legend.position = "bottom", 
            axis.line = element_blank(), 
            legend.justification = "center", 
            plot.title = element_text(hjust = 0.5), 
            axis.title.y = element_text(angle = 360, vjust = 0.5), 
            panel.border = element_rect(fill = NA, size = 1, color = "black"), 
            axis.text.x = element_text(angle = 45, size = 8, vjust = 1, hjust = 1)) + 
      guides(color = guide_colorbar(title.position = "top", title.hjust = 0.5, barwidth = unit(5, "cm"), order = 2), 
             size = guide_legend(title.position = "top", title.hjust = 0.5, order = 1))

We can see C1QA, which is a myeloid marker, in cluster 3. Cluster 4 shows CD3D & NKG7; this indicates that the cluster is composed of T/NK cells. Cluster 2 has TOP2A, which is a DC marker.

p5

Cell Type Identification

SingleR: scRNA-seq Reference

SingleR uses Spearman correlation to compare scRNA-seq data to a labeled reference dataset. We’ll use the annotated cells from Baron et al (2016) in order to give broad cell types to our clusters. Note: these should not be considered ground truth, but can help guide our manual annotations. Also, this dataset can be loaded from the scRNAseq package, but I’ve created a version normalized using SCTransform instead og the default log-normalization, to better match the cells we’d like to annotate.

Primary Tumor

We read in the SCTransform-normalized reference dataset, and run SingleR.

sc_ref <- readRDS("/Volumes/labs/Home/Jen Jen Yeh Lab/Jack/scRNAseq/Seurat/single_cell_ref_normalized.Rds")
cell_preds <- SingleR(test = primary@assays$integrated@scale.data, 
                      ref = sc_ref, 
                      labels = sc_ref$label, 
                      clusters = primary$seurat_clusters, 
                      de.method = "wilcox")
primary[["SingleR_sc"]] <- cell_preds$labels[match(primary[[]][["seurat_clusters"]], rownames(cell_preds))]
primary[["SingleR_sc"]] <- case_when(primary[["SingleR_sc"]] == "acinar" ~ "Acinar", 
                                     primary[["SingleR_sc"]] == "activated_stellate" ~ "Activated Stellate", 
                                     primary[["SingleR_sc"]] == "ductal" ~ "Epithelial", 
                                     primary[["SingleR_sc"]] == "endothelial" ~ "Endothelial", 
                                     primary[["SingleR_sc"]] == "macrophage" ~ "Macrophage", 
                                     primary[["SingleR_sc"]] == "mast" ~ "Mast", 
                                     primary[["SingleR_sc"]] == "quiescent_stellate" ~ "Quiescent Stellate")
p6 <- DimPlot(primary, group.by = "SingleR_sc") + 
      labs(x = "UMAP 1", y = "UMAP 2", subtitle = "Primary Tumor") + 
      scale_color_paletteer_d("ggthemes::Classic_20") + 
      theme_yehlab() + 
      theme(legend.position = "right", 
            plot.title = element_blank(), 
            legend.text = element_text(size = 10), 
            plot.subtitle = element_text(), 
            axis.title = element_text(hjust = 0, size = 9)) + 
      guides(color = guide_legend(ncol = 1, override.aes = list(size = 3)))

Metastatic Tumor

We’ll do the same for the metastatic tumor cells.

cell_preds <- SingleR(test = meta@assays$integrated@scale.data, 
                      ref = sc_ref, 
                      labels = sc_ref$label, 
                      clusters = meta$seurat_clusters, 
                      de.method = "wilcox")
meta[["SingleR_sc"]] <- cell_preds$labels[match(meta[[]][["seurat_clusters"]], rownames(cell_preds))]
meta[["SingleR_sc"]] <- case_when(meta[["SingleR_sc"]] == "acinar" ~ "Acinar", 
                                  meta[["SingleR_sc"]] == "t_cell" ~ "T", 
                                  meta[["SingleR_sc"]] == "endothelial" ~ "Endothelial", 
                                  meta[["SingleR_sc"]] == "macrophage" ~ "Macrophage", 
                                  meta[["SingleR_sc"]] == "gamma" ~ "Gamma")
p7 <- DimPlot(meta, group.by = "SingleR_sc", cols = paletteer_d("ggthemes::Classic_20")[8:12]) + 
      labs(subtitle = "Metastatic Tumor") + 
      theme_yehlab() + 
      theme(legend.position = "right", 
            plot.title = element_blank(), 
            axis.title = element_blank(), 
            legend.text = element_text(size = 9), 
            plot.subtitle = element_text()) + 
      guides(color = guide_legend(ncol = 1, override.aes = list(size = 3)))

Let’s check out the results! I really doubt the Gamma label for the metastatic samples, as Gamma cells make up around 3-5% of all pancreatic islet cells. The T cells and Macrophage labels match up with the authors’ annotations. I would guess that the rest of the cells are split between normal ductal cells and PDAC cells.

As for the primary tumor cells, we see a Macrophage cluster and an Activated Stellate (read: CAF) cluster which matches the authors’ annotations. The acinar cells are interesting, as the authors noted that one of the tumor cluster was enriched for acinar marker genes. Most likely though, they are truly Epithelial cells which we’ll investigate further. Finally, the small Endothelial cluster roughly matches the size of the Endothelial cluster identified by the authors, though here it’s located next to the putative Epithelial cells instead of next to the CAFs as it was in the original publication.

(p7 / p6) + plot_annotation(title = "SingleR Annotations", subtitle = "scRNA-seq Reference", 
                            theme = theme(plot.title = element_text(size = 16, hjust = 0.5), 
                                          plot.subtitle = element_text(face = "italic", hjust = 0.5, size = 9)))

SingleR: Bulk RNA-seq Reference

We’ll perform the same procedure, this time using a bulk RNA-seq reference - data from the Human Primary Cell Atlas.

bulk_ref <- HumanPrimaryCellAtlasData()

Primary Tumor

We run SingleR again, this time using the bulk reference.

cell_preds <- SingleR(test = primary@assays$integrated@scale.data, 
                      ref = bulk_ref, 
                      labels = bulk_ref$label.main, 
                      clusters = primary$seurat_clusters, 
                      de.method = "wilcox")
primary[["SingleR_bulk"]] <- cell_preds$labels[match(primary[[]][["seurat_clusters"]], rownames(cell_preds))]
primary[["SingleR_bulk"]] <- case_when(primary[["SingleR_bulk"]] == "Epithelial_cells" ~ "Epithelial", 
                                       primary[["SingleR_bulk"]] == "Endothelial_cells" ~ "Endothelial", 
                                       primary[["SingleR_bulk"]] == "Fibroblasts" ~ "Fibroblast", 
                                       primary[["SingleR_bulk"]] == "Neutrophils" ~ "Neutrophil", 
                                       primary[["SingleR_bulk"]] == "Macrophage" ~ "Macrophage", 
                                       primary[["SingleR_bulk"]] == "BM & Prog." ~ "Bone Marrow Progenitor")
p8 <- DimPlot(primary, group.by = "SingleR_bulk") + 
      labs(x = "UMAP 1", y = "UMAP 2", subtitle = "Primary Tumor") + 
      scale_color_paletteer_d("ggthemes::Classic_20") + 
      theme_yehlab() + 
      theme(legend.position = "right", 
            plot.title = element_blank(), 
            legend.text = element_text(size = 9), 
            plot.subtitle = element_text(), 
            axis.title = element_text(hjust = 0, size = 9)) + 
      guides(color = guide_legend(ncol = 1, override.aes = list(size = 3)))

Metastatic Tumor

We’ll do the same for the metastatic tumor cells.

cell_preds <- SingleR(test = meta@assays$integrated@scale.data, 
                      ref = bulk_ref, 
                      labels = bulk_ref$label.main, 
                      clusters = meta$seurat_clusters, 
                      de.method = "wilcox")
meta[["SingleR_bulk"]] <- cell_preds$labels[match(meta[[]][["seurat_clusters"]], rownames(cell_preds))]
meta[["SingleR_bulk"]] <- case_when(meta[["SingleR_bulk"]] == "Epithelial_cells" ~ "Epithelial", 
                                  meta[["SingleR_bulk"]] == "Monocyte" ~ "Monocyte", 
                                  meta[["SingleR_bulk"]] == "T_cells" ~ "T", 
                                  meta[["SingleR_bulk"]] == "Neurons" ~ "Neuron")
p9 <- DimPlot(meta, group.by = "SingleR_bulk", cols = paletteer_d("ggthemes::Classic_20")[8:11]) + 
      labs(subtitle = "Metastatic Tumor") + 
      theme_yehlab() + 
      theme(legend.position = "right", 
            plot.title = element_blank(), 
            axis.title = element_blank(), 
            legend.text = element_text(size = 9), 
            plot.subtitle = element_text()) + 
      guides(color = guide_legend(ncol = 1, override.aes = list(size = 3)))

Let’s check out the results! For the primary tumor sample: we can be fairly certain that the large Epithelial cluster identity is correct, as we saw it in the single cell and bulk reference annotations. I think the Bone Marrow Progenitor cluster is most likely the Epithelial \(\rightarrow\) Mesenchymal cluster the authors annotated. The Fibroblast & Endothelial cluster identities also line up with what we saw in the single cell annotations. Lastly, I’m interested in the annotation of the Neutrophil cluster, as that’s something the authors did not annotate.

For the metastatic samples, we see again the small T & Macrophage clusters. We can be pretty sure that those are correct. I trust that the large Epithelial cluster is correct, though I’m puzzled as to why the final cluster was annotated as containing Neurons. The single cell reference defined that cluster as Endothelial. We’ll have to do some further investigation on that group of cells.

(p9 / p8) + plot_annotation(title = "SingleR Annotations", subtitle = "Bulk RNA-seq Reference", 
                            theme = theme(plot.title = element_text(size = 16, hjust = 0.5), 
                                          plot.subtitle = element_text(face = "italic", hjust = 0.5, size = 9)))

Pathway Analysis

Primary

We’ll use Dr. Xianlu Peng’s DECODER method in order to perform a pancreas-specific NMF version of pathway analysis, in which each cell is assigned compartment weights for pathways such as Immune, Classical PDAC, etc. The PDAC subtypes used here, Basal & Classical, are taken from Moffitt et al (2015).

sample_wts_unscaled <- Decon_single_sample("TCGA_RNAseq_PAAD", 
                                           primary@assays$integrated@data, 
                                           "geneSymbol")
sample_wts <- Norm_PDAC_weights(sample_wts_unscaled)
primary <- AddMetaData(primary, sample_wts$Immune, "immune")
primary <- AddMetaData(primary, sample_wts$bcRatio, "bc_ratio")
primary <- AddMetaData(primary, sample_wts$Endocrine, "endocrine")
primary <- AddMetaData(primary, sample_wts_unscaled[, 9], "basal")
primary <- AddMetaData(primary, sample_wts_unscaled[, 5], "classical")
primary <- AddMetaData(primary, sample_wts$NormalStroma, "norm_stroma")
primary <- AddMetaData(primary, sample_wts$ActivatedStroma, "act_stroma")

Now let’s generate some plots! We see that the Immune & Fibroblast scores correspond with what we saw in SingleR. Interestingly, it appears as though the PDAC scores are high in only one of the Epithelial clusters (Cluster 1). This leads me to believe that there are potential subclusters within the cluster, and that the other cluster is composed of normal, non-malignant Epithelial cells.

plots <- list()
vars <- c("immune", "endocrine", "basal", "classical", "norm_stroma", "act_stroma")
titles <- c("Immune", "Endocrine", "Basal PDAC", "Classical PDAC", "Normal Stroma", "Activated Stroma")
for (i in seq_along(vars)) {
  plots[[i]] <- FeaturePlot(primary, features = vars[i], cols = c("grey70", "firebrick")) + 
                labs(subtitle = titles[i]) + 
                theme_yehlab() + 
                theme(legend.position = "none", 
                      axis.title = element_blank(), 
                      plot.title = element_blank(), 
                      plot.title.position = "plot", 
                      plot.subtitle = element_text(hjust = 0.5, size = 12))
}
p10 <- ggarrange(plotlist = plots, ncol = 3, nrow = 2) %>% 
       annotate_figure(bottom = "UMAP 1", left = text_grob("UMAP 2", rot = 360), 
                       top = "DECODER Pathway Analysis (Primary Tumor)")
p10

Metastatic Tumor

We’ll repeat the DECODER analysis and plot generation for the metastatic cells.

sample_wts_unscaled <- Decon_single_sample("TCGA_RNAseq_PAAD", 
                                           meta@assays$integrated@data, 
                                           "geneSymbol")
sample_wts <- Norm_PDAC_weights(sample_wts_unscaled)
meta <- AddMetaData(meta, sample_wts$Immune, "immune")
meta <- AddMetaData(meta, sample_wts$bcRatio, "bc_ratio")
meta <- AddMetaData(meta, sample_wts$Endocrine, "endocrine")
meta <- AddMetaData(meta, sample_wts_unscaled[, 9], "basal")
meta <- AddMetaData(meta, sample_wts_unscaled[, 5], "classical")
meta <- AddMetaData(meta, sample_wts$NormalStroma, "norm_stroma")
meta <- AddMetaData(meta, sample_wts$ActivatedStroma, "act_stroma")

Now let’s generate some plots!

plots <- list()
vars <- c("immune", "endocrine", "basal", "classical", "norm_stroma", "act_stroma")
titles <- c("Immune", "Endocrine", "Basal PDAC", "Classical PDAC", "Normal Stroma", "Activated Stroma")
for (i in seq_along(vars)) {
  plots[[i]] <- FeaturePlot(meta, features = vars[i], cols = c("grey70", "firebrick")) + 
                labs(subtitle = titles[i]) + 
                theme_yehlab() + 
                theme(legend.position = "none", 
                      axis.title = element_blank(), 
                      plot.title = element_blank(), 
                      plot.title.position = "plot", 
                      plot.subtitle = element_text(hjust = 0.5, size = 12))
}

Despite some missing genes, we see that our T & Macrophage annotations match the Immune weights, and that the Basal & Classical PDAC scores are located in the putative Epithelial clusters. Specifically, the Basal PDAC weights are highest in the odd cluster that was annotated as Gamma / Neuron cells by SingleR. All three Epithelial clusters, though, do have either Basal or Classical PDAC cells in them.

p11 <- ggarrange(plotlist = plots, ncol = 3, nrow = 2) %>% 
       annotate_figure(bottom = "UMAP 1", left = text_grob("UMAP 2", rot = 360), 
                       top = "DECODER Pathway Analysis (Metastatic Tumor)")
p11

Assign Broad Cell Types

We have enough information to accurately assign broad cell types to our clusters.

Primary Tumor

primary$cell_type <- case_when(primary$seurat_clusters == 0 ~ "Epithelial", 
                               primary$seurat_clusters == 1 ~ "Epithelial", 
                               primary$seurat_clusters == 2 ~ "Fibroblast", 
                               primary$seurat_clusters == 3 ~ "Macrophage", 
                               primary$seurat_clusters == 4 ~ "Neutrophil", 
                               primary$seurat_clusters == 5 ~ "Mesenchymal Stem Cell", 
                               primary$seurat_clusters == 6 ~ "Fibroblast", 
                               primary$seurat_clusters == 7 ~ "Endothelial")

Metastatic Tumor

meta$cell_type <- case_when(meta$seurat_clusters == 0 ~ "Epithelial", 
                            meta$seurat_clusters == 1 ~ "Epithelial", 
                            meta$seurat_clusters == 2 ~ "Epithelial", 
                            meta$seurat_clusters == 3 ~ "Macrophage", 
                            meta$seurat_clusters == 4 ~ "T")

CNV Estimation

Next we’ll attempt to identify malignant cells using single-cell copy number variation estimation as implemented in the CONCISmat package. Details of the GMM methodology used can be found at the Diaz Lab’s GitHub repository.

Primary Tumor

chrom_regions <- read.table("/Volumes/labs/Home/Jen Jen Yeh Lab/Jack/scRNAseq/chrom_arm_positions.txt", 
                            sep = "\t", 
                            row.names = 1, 
                            header = TRUE)
gene_pos <- getGenePositions(rownames(primary))
cpm <- t(t(as.matrix(primary@assays$SCT@counts)) / colSums(as.matrix(primary@assays$SCT@counts))) * 10^5
cpm <- log2(cpm + 1)
norm_factor <- calcNormFactors(cpm)
cnv_est <- plotAll(mat = cpm, 
                   normFactor = norm_factor, 
                   regions = chrom_regions, 
                   gene_pos = gene_pos, 
                   fname = "../Data/Primary")

## [1] "Fitting GMM for chr1 0:122026459 iteration 1"
## number of iterations= 512 
## [1] "Fitting GMM for chr1 0:122026459 iteration 2"
## number of iterations= 405 
## [1] "Fitting GMM for chr1 0:122026459 iteration 3"
## number of iterations= 354 
## [1] "Fitting GMM for chr1 0:122026459 iteration 4"
## number of iterations= 499 
## [1] 44

## [1] "Fitting GMM for chr1 124932724:248956422 iteration 1"
## number of iterations= 438 
## [1] "Fitting GMM for chr1 124932724:248956422 iteration 2"
## number of iterations= 236 
## [1] "Fitting GMM for chr1 124932724:248956422 iteration 3"
## number of iterations= 331 
## [1] "Fitting GMM for chr1 124932724:248956422 iteration 4"
## number of iterations= 221 
## [1] 44

## [1] "Fitting GMM for chr2 0:92188145 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr2 0:92188145 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr2 0:92188145 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr2 0:92188145 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr2 94090557:242193529 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr2 94090557:242193529 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr2 94090557:242193529 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr2 94090557:242193529 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr3 0:90772458 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr3 0:90772458 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr3 0:90772458 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr3 0:90772458 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr3 93655574:198295559 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr3 93655574:198295559 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr3 93655574:198295559 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr3 93655574:198295559 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr4 51743951:190214555 iteration 1"
## number of iterations= 634 
## [1] "Fitting GMM for chr4 51743951:190214555 iteration 2"
## number of iterations= 371 
## [1] "Fitting GMM for chr4 51743951:190214555 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr4 51743951:190214555 iteration 4"
## number of iterations= 943 
## [1] 44

## [1] "Fitting GMM for chr5 50059807:181538259 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr5 50059807:181538259 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr5 50059807:181538259 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr5 50059807:181538259 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr7 61528020:159345973 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr7 61528020:159345973 iteration 2"
## number of iterations= 947 
## [1] "Fitting GMM for chr7 61528020:159345973 iteration 3"
## number of iterations= 669 
## [1] "Fitting GMM for chr7 61528020:159345973 iteration 4"
## number of iterations= 573 
## [1] 44

## [1] "Fitting GMM for chr9 45518558:138394717 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr9 45518558:138394717 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr9 45518558:138394717 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr9 45518558:138394717 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr10 41593521:133797422 iteration 1"
## number of iterations= 761 
## [1] "Fitting GMM for chr10 41593521:133797422 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr10 41593521:133797422 iteration 3"
## number of iterations= 509 
## [1] "Fitting GMM for chr10 41593521:133797422 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr11 54425074:135086622 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr11 54425074:135086622 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr11 54425074:135086622 iteration 3"
## number of iterations= 848 
## [1] "Fitting GMM for chr11 54425074:135086622 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr12 37185252:133275309 iteration 1"
## number of iterations= 473 
## [1] "Fitting GMM for chr12 37185252:133275309 iteration 2"
## number of iterations= 674 
## [1] "Fitting GMM for chr12 37185252:133275309 iteration 3"
## number of iterations= 668 
## [1] "Fitting GMM for chr12 37185252:133275309 iteration 4"
## number of iterations= 561 
## [1] 44

## [1] "Fitting GMM for chr14 18173523:107043718 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr14 18173523:107043718 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr14 18173523:107043718 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr14 18173523:107043718 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr15 19725254:101991189 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr15 19725254:101991189 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr15 19725254:101991189 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr15 19725254:101991189 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr17 26566633:83257441 iteration 1"
## number of iterations= 675 
## [1] "Fitting GMM for chr17 26566633:83257441 iteration 2"
## number of iterations= 654 
## [1] "Fitting GMM for chr17 26566633:83257441 iteration 3"
## number of iterations= 538 
## [1] "Fitting GMM for chr17 26566633:83257441 iteration 4"
## number of iterations= 684 
## [1] 44

## [1] "Fitting GMM for chr19 0:24498980 iteration 1"
## number of iterations= 562 
## [1] "Fitting GMM for chr19 0:24498980 iteration 2"
## number of iterations= 633 
## [1] "Fitting GMM for chr19 0:24498980 iteration 3"
## number of iterations= 661 
## [1] "Fitting GMM for chr19 0:24498980 iteration 4"
## number of iterations= 505 
## [1] 44

## [1] "Fitting GMM for chr19 27190874:58617616 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr19 27190874:58617616 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr19 27190874:58617616 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr19 27190874:58617616 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

bic_table <- read.table("../Data/Primary_BIC_LR.txt", 
                        sep = "\t", 
                        row.names = 1, 
                        header = TRUE, 
                        check.names = FALSE)
cand_regions <- rownames(bic_table[bic_table$`BIC difference` > 200 & bic_table$`LRT adj. p-val` < .01, ])

We use \(k = 3\) clusters, as we assume that the cells can be divided into 1) PDAC cells, 2) CAF cells, and 3) non-malignant cells. We see that chromosomes 9q, 12q, and 17q have the most estimated CNVs.

hist1 <- plotHistogram(cnv_est[, cand_regions], 
                       cpm, 
                       clusters = 3, 
                       zscoreThreshold = 3, 
                       celltypes = primary$cell_type, 
                       patients = primary$sample)

normal_primary <- which(hist1 == 1)
primary@meta.data$CNV <- ifelse(rownames(primary@meta.data) %in% names(normal_primary), 
                                "Normal", "Malignant")
p12 <- DimPlot(primary, group.by = "CNV") + 
       labs(subtitle = "Primary Tumor") + 
       theme_yehlab() + 
       theme(legend.position = "none", 
             plot.title = element_blank(), 
             axis.title = element_blank(), 
             plot.subtitle = element_text(hjust = 0.5))

Metastatic Tumor

We repeat the process for the metastatic tumor samples.

gene_pos <- getGenePositions(rownames(meta))
cpm <- t(t(as.matrix(meta@assays$SCT@counts)) / colSums(as.matrix(meta@assays$SCT@counts))) * 10^5
cpm <- log2(cpm + 1)
norm_factor <- calcNormFactors(cpm)
cnv_est <- plotAll(mat = cpm, 
                   normFactor = norm_factor, 
                   regions = chrom_regions, 
                   gene_pos = gene_pos, 
                   fname = "../Data/Meta")

## [1] "Fitting GMM for chr1 0:122026459 iteration 1"
## number of iterations= 526 
## [1] "Fitting GMM for chr1 0:122026459 iteration 2"
## number of iterations= 504 
## [1] "Fitting GMM for chr1 0:122026459 iteration 3"
## number of iterations= 487 
## [1] "Fitting GMM for chr1 0:122026459 iteration 4"
## number of iterations= 283 
## [1] 44

## [1] "Fitting GMM for chr1 124932724:248956422 iteration 1"
## number of iterations= 413 
## [1] "Fitting GMM for chr1 124932724:248956422 iteration 2"
## number of iterations= 563 
## [1] "Fitting GMM for chr1 124932724:248956422 iteration 3"
## number of iterations= 304 
## [1] "Fitting GMM for chr1 124932724:248956422 iteration 4"
## number of iterations= 538 
## [1] 44

## [1] "Fitting GMM for chr2 0:92188145 iteration 1"
## number of iterations= 807 
## [1] "Fitting GMM for chr2 0:92188145 iteration 2"
## number of iterations= 660 
## [1] "Fitting GMM for chr2 0:92188145 iteration 3"
## number of iterations= 839 
## [1] "Fitting GMM for chr2 0:92188145 iteration 4"
## number of iterations= 651 
## [1] 44

## [1] "Fitting GMM for chr2 94090557:242193529 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr2 94090557:242193529 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr2 94090557:242193529 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr2 94090557:242193529 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr3 93655574:198295559 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr3 93655574:198295559 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr3 93655574:198295559 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr3 93655574:198295559 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr4 51743951:190214555 iteration 1"
## number of iterations= 959 
## [1] "Fitting GMM for chr4 51743951:190214555 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr4 51743951:190214555 iteration 3"
## number of iterations= 997 
## [1] "Fitting GMM for chr4 51743951:190214555 iteration 4"
## number of iterations= 718 
## [1] 44

## [1] "Fitting GMM for chr5 50059807:181538259 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr5 50059807:181538259 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr5 50059807:181538259 iteration 3"
## number of iterations= 837 
## [1] "Fitting GMM for chr5 50059807:181538259 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr6 0:58553888 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr6 0:58553888 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr6 0:58553888 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr6 0:58553888 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr7 61528020:159345973 iteration 1"
## number of iterations= 657 
## [1] "Fitting GMM for chr7 61528020:159345973 iteration 2"
## number of iterations= 359 
## [1] "Fitting GMM for chr7 61528020:159345973 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr7 61528020:159345973 iteration 4"
## number of iterations= 715 
## [1] 44

## [1] "Fitting GMM for chr9 45518558:138394717 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr9 45518558:138394717 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr9 45518558:138394717 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr9 45518558:138394717 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

## [1] "Fitting GMM for chr10 41593521:133797422 iteration 1"
## number of iterations= 526 
## [1] "Fitting GMM for chr10 41593521:133797422 iteration 2"
## number of iterations= 274 
## [1] "Fitting GMM for chr10 41593521:133797422 iteration 3"
## number of iterations= 557 
## [1] "Fitting GMM for chr10 41593521:133797422 iteration 4"
## number of iterations= 372 
## [1] 44

## [1] "Fitting GMM for chr11 54425074:135086622 iteration 1"
## number of iterations= 651 
## [1] "Fitting GMM for chr11 54425074:135086622 iteration 2"
## number of iterations= 441 
## [1] "Fitting GMM for chr11 54425074:135086622 iteration 3"
## number of iterations= 788 
## [1] "Fitting GMM for chr11 54425074:135086622 iteration 4"
## number of iterations= 487 
## [1] 44

## [1] "Fitting GMM for chr12 37185252:133275309 iteration 1"
## number of iterations= 133 
## [1] "Fitting GMM for chr12 37185252:133275309 iteration 2"
## number of iterations= 275 
## [1] "Fitting GMM for chr12 37185252:133275309 iteration 3"
## number of iterations= 148 
## [1] "Fitting GMM for chr12 37185252:133275309 iteration 4"
## number of iterations= 286 
## [1] 44

## [1] "Fitting GMM for chr14 18173523:107043718 iteration 1"
## number of iterations= 788 
## [1] "Fitting GMM for chr14 18173523:107043718 iteration 2"
## number of iterations= 385 
## [1] "Fitting GMM for chr14 18173523:107043718 iteration 3"
## number of iterations= 531 
## [1] "Fitting GMM for chr14 18173523:107043718 iteration 4"
## number of iterations= 456 
## [1] 44

## [1] "Fitting GMM for chr15 19725254:101991189 iteration 1"
## number of iterations= 485 
## [1] "Fitting GMM for chr15 19725254:101991189 iteration 2"
## number of iterations= 351 
## [1] "Fitting GMM for chr15 19725254:101991189 iteration 3"
## number of iterations= 267 
## [1] "Fitting GMM for chr15 19725254:101991189 iteration 4"
## number of iterations= 322 
## [1] 44

## [1] "Fitting GMM for chr17 26566633:83257441 iteration 1"
## number of iterations= 350 
## [1] "Fitting GMM for chr17 26566633:83257441 iteration 2"
## number of iterations= 282 
## [1] "Fitting GMM for chr17 26566633:83257441 iteration 3"
## number of iterations= 379 
## [1] "Fitting GMM for chr17 26566633:83257441 iteration 4"
## number of iterations= 341 
## [1] 44

## [1] "Fitting GMM for chr19 0:24498980 iteration 1"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr19 0:24498980 iteration 2"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr19 0:24498980 iteration 3"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] "Fitting GMM for chr19 0:24498980 iteration 4"
## number of iterations= 220 
## [1] 44

## [1] "Fitting GMM for chr19 27190874:58617616 iteration 1"
## number of iterations= 423 
## [1] "Fitting GMM for chr19 27190874:58617616 iteration 2"
## number of iterations= 429 
## [1] "Fitting GMM for chr19 27190874:58617616 iteration 3"
## number of iterations= 470 
## [1] "Fitting GMM for chr19 27190874:58617616 iteration 4"
## WARNING! NOT CONVERGENT! 
## number of iterations= 1000 
## [1] 44

bic_table <- read.table("../Data/Meta_BIC_LR.txt", 
                        sep = "\t", 
                        row.names = 1, 
                        header = TRUE, 
                        check.names = FALSE)
cand_regions <- rownames(bic_table[bic_table$`BIC difference` > 200 & bic_table$`LRT adj. p-val` < .01, ])

We see that chromosomes 2p, 7q, 12q, and 15q have the most estimated CNVs for the metastatic samples.

hist2 <- plotHistogram(cnv_est[, cand_regions], 
                       cpm, 
                       clusters = 3, 
                       zscoreThreshold = 3, 
                       celltypes = meta$cell_type, 
                       patients = meta$sample)

normal_meta <- which(hist2 == 3)
meta@meta.data$CNV <- ifelse(rownames(meta@meta.data) %in% names(normal_meta), 
                             "Normal", "Malignant")
p13 <- DimPlot(meta, group.by = "CNV") + 
       labs(subtitle = "Metastatic Tumor") + 
       theme_yehlab() + 
       theme(legend.position = "none", 
             plot.title = element_blank(), 
             axis.title = element_blank(), 
             plot.subtitle = element_text(hjust = 0.5))

Let’s check out the results! We see that the cells estimated to be malignant are mostly located in the Epithelial clusters. This is what we’d expect from Pancreatic Ductal Adenocarcinoma. There are also a decent amount of malignant cells in the MSC cluster we identified. Notably, we don’t see a significant concentration of malignant cells in the Fibroblast cluster from the primary samples. This makes me doubt the authors’ annotation of those cells as Cancer Associated Fibroblasts (CAFs).

ggarrange(p12, p13, ncol = 2) %>% 
  annotate_figure(bottom = "UMAP 1", left = text_grob("UMAP 2", rot = 360), 
                  top = "Malignant Cells Estimated Through CONICSmat")

Reclustering Cells

Primary Tumor - Epithelial Cells

Lastly, we’ll use my method, SCISSORS, to determine whether separate Basal & Classical PDAC clusters exist in the PDAC cluster (Cluster 1) present in the primary tumor samples. The method re-processes the data and tests several clustering parameter sets in order to determine which parameter set fits the cells the best. This is done by maximizing the mean silhouette score for each cluster. Here we run the method and generate some plots of the results.

epi_reclust <- ReclusterCells(seurat.object = primary, 
                              which.clust = 1, 
                              n.HVG = 3000, 
                              n.PC = 15, 
                              k.vals = c(20, 30, 40), 
                              resolution.vals = c(.2, .3, .4), 
                              redo.embedding = TRUE, 
                              random.seed = 629)

## [1] "Reclustering cells in cluster 1 using k = 40 & resolution = 0.2; S = 0.349"

p14 <- DimPlot(epi_reclust) + 
       labs(subtitle = "SCISSORS - Cluster 1") + 
       theme_yehlab() + 
       theme(legend.position = "none", 
             plot.subtitle = element_text(hjust = 0.5), 
             axis.title = element_blank()) + 
       guides(color = guide_legend(ncol = 1, override.aes = list(size = 3)))
p15 <- FeaturePlot(epi_reclust, features = "classical", cols = c("grey90", "firebrick")) + 
       labs(subtitle = "Classical PDAC") + 
       theme_yehlab() + 
       theme(legend.position = "none",
             axis.title = element_blank(), 
             plot.title = element_blank(), 
             plot.subtitle = element_text(hjust = 0.5))
p16 <- FeaturePlot(epi_reclust, features = "basal", cols = c("grey90", "firebrick")) + 
       labs(subtitle = "Basal PDAC") + 
       theme_yehlab() + 
       theme(legend.position = "none",
             axis.title = element_blank(), 
             plot.title = element_blank(), 
             plot.subtitle = element_text(hjust = 0.5))
p0a <- p0 + 
       labs(subtitle = "Primary Tumor") + 
       theme(axis.title.y = element_blank(), 
             axis.title.x = element_blank(), 
             plot.title = element_blank(), 
             legend.position = "right", 
             plot.subtitle = element_text(hjust = 0.5)) + 
       guides(color = guide_legend(ncol = 1, override.aes = list(size = 3)))
p17 <- (p14 | p0a) / (p15 | p16)

We see that while it’s possible to find subclusters within Cluster 1, the “best” clustering does not separate out the Basal & Classical components of the PDAC cluster.

ggarrange(p17) %>% 
  annotate_figure(bottom = "UMAP 1", left = text_grob("UMAP 2", rot = 360))

Primary Tumor - Fibroblasts

Lastly, we saw in the DECODER results that one corner of the Fibroblast cluster in the primary tumor cells had high Endocrine weights. We’ll run SCISSORS on that cluster (Cluster 2) in order to see if subgroups can be extracted.

fibro_reclust <- ReclusterCells(seurat.object = primary, 
                                which.clust = 2, 
                                n.HVG = 3000, 
                                n.PC = 15, 
                                k.vals = c(10, 20, 30), 
                                resolution.vals = c(.2, .3, .4), 
                                redo.embedding = TRUE, 
                                random.seed = 629)

## [1] "Reclustering cells in cluster 2 using k = 30 & resolution = 0.2; S = 0.441"

fibro_markers <- FindAllMarkers(fibro_reclust, 
                                test.use = "wilcox", 
                                logfc.threshold = 2, 
                                only.pos = TRUE, 
                                verbose = FALSE, 
                                random.seed = 629)

We generate some plots of the results.

p18 <- DimPlot(fibro_reclust) + 
       labs(subtitle = "SCISSORS - Cluster 2") + 
       theme_yehlab() + 
       theme(legend.position = "none", 
             plot.subtitle = element_text(hjust = 0.5), 
             axis.title = element_blank()) + 
       guides(color = guide_legend(ncol = 1, override.aes = list(size = 3)))
p19 <- FeaturePlot(fibro_reclust, features = "act_stroma", cols = c("grey90", "firebrick")) + 
       labs(subtitle = "Activated Stroma") + 
       theme_yehlab() + 
       theme(legend.position = "none",
             axis.title = element_blank(), 
             plot.title = element_blank(), 
             plot.subtitle = element_text(hjust = 0.5))
p20 <- FeaturePlot(fibro_reclust, features = "endocrine", cols = c("grey90", "firebrick")) + 
       labs(subtitle = "Endocrine") + 
       theme_yehlab() + 
       theme(legend.position = "none",
             axis.title = element_blank(), 
             plot.title = element_blank(), 
             plot.subtitle = element_text(hjust = 0.5))
p21 <- (p18 | p0a) / (p19 | p20)

We see two clusters - one of which clearly corresponds to Activated Stroma / CAF cells, and the other of which seems so be composed of islet cells, as they have high endocrine pancreas scores.

ggarrange(p21) %>% 
  annotate_figure(bottom = "UMAP 1", left = text_grob("UMAP 2", rot = 360))

The top 5 markers for each cluster are a mixture of proto-oncogenes and stroma-related genes.

fibro_markers %>% 
  dplyr::select(cluster, gene, avg_log2FC, p_val_adj, pct.1, pct.2) %>% 
  group_by(cluster) %>% 
  slice_head(n = 5) %>% 
  kbl(booktabs = TRUE, col.names = c("Cluster", "Gene", "Mean Log2FC", 
                                     "Adj. P-value", "% Expressed 1", "% Expressed 2")) %>% 
  kable_minimal("hover", full_width = FALSE)

Cluster	Gene	Mean Log2FC	Adj. P-value	% Expressed 1	% Expressed 2
0	ITGB1	3.955181	0	0.963	0.564
0	TPM4	2.875344	0	0.961	0.569
0	SDC1	3.867239	0	0.809	0.157
0	MYL6	3.575856	0	0.972	0.795
0	IFI27	3.880991	0	0.770	0.182
1	JUNB	4.690781	0	0.986	0.594
1	JUN	5.476523	0	0.930	0.404
1	IER2	3.551652	0	0.977	0.504
1	FOSB	3.569275	0	0.861	0.287
1	FOS	3.576318	0	0.986	0.687

Assign Cell Types

Finally, we have enough information to accurately assign cell types to our clusters.

Primary Tumor

endocrine <- names(fibro_reclust$seurat_clusters[fibro_reclust$seurat_clusters == 1])
CAF <- names(fibro_reclust$seurat_clusters[fibro_reclust$seurat_clusters == 0])
primary$cell_type <- case_when(rownames(primary@meta.data) %in% endocrine ~ "Endocrine", 
                               rownames(primary@meta.data) %in% CAF ~ "CAF", 
                               primary$seurat_clusters == 1 ~ "PDAC", 
                               TRUE ~ primary$cell_type)
Idents(primary) <- "cell_type"
p22 <- DimPlot(primary, group.by = "cell_type") + 
       scale_color_paletteer_d("ggthemes::Classic_20") + 
       labs(subtitle = "Primary Tumor", x = "UMAP 1", y = "UMAP 2") + 
       theme_yehlab() + 
       theme(legend.position = "right", 
             plot.title = element_blank(), 
             legend.text = element_text(size = 10), 
             plot.subtitle = element_text(), 
             axis.title = element_text(hjust = 0, size = 9)) + 
       guides(color = guide_legend(ncol = 1, override.aes = list(size = 3)))

Metastatic Tumor

meta$cell_type <- case_when(meta$seurat_clusters == 0 ~ "PDAC", 
                            meta$seurat_clusters == 1 ~ "PDAC", 
                            meta$seurat_clusters == 2 ~ "PDAC", 
                            meta$seurat_clusters == 3 ~ "Macrophage", 
                            meta$seurat_clusters == 4 ~ "T")
Idents(meta) <- "cell_type"
p23 <- DimPlot(meta, group.by = "cell_type", cols = paletteer_d("ggthemes::Classic_20")[c(6, 9, 11)]) + 
       labs(subtitle = "Metastatic Tumor") + 
       theme_yehlab() + 
       theme(legend.position = "right", 
             plot.title = element_blank(), 
             axis.title = element_blank(), 
             legend.text = element_text(size = 9), 
             plot.subtitle = element_text()) + 
       guides(color = guide_legend(ncol = 1, override.aes = list(size = 3)))

Let’s check out the final results!

(p23 / p22) + plot_annotation(title = "Final Celltype Labels", 
                              theme = theme(plot.title = element_text(size = 16, hjust = 0.5)))

Differential Expression

Lastly, we’ll run differential expression tests again using our final cell labels this time instead of our cluster identities.

Primary Tumor

primary_markers_final <- FindAllMarkers(primary, 
                                        logfc.threshold = 2, 
                                        test.use = "wilcox", 
                                        verbose = FALSE, 
                                        only.pos = TRUE, 
                                        random.seed = 629)
primary_markers_final_top5 <- primary_markers_final %>% 
                              group_by(cluster) %>% 
                              slice_head(n = 5)
p24 <- DotPlot(primary, features = primary_markers_final_top5$gene, cols = c("grey90", "firebrick")) + 
       labs(x = "Gene", 
            y = "Cell", 
            color = "Mean Expression", 
            size = "Percent Expressed", 
            title = "Primary Tumor Marker Genes") + 
       theme(legend.position = "bottom", 
             axis.line = element_blank(), 
             legend.justification = "center", 
             plot.title = element_text(hjust = 0.5), 
             axis.title.y = element_blank(), 
             panel.border = element_rect(fill = NA, size = 1, color = "black"), 
             axis.text.x = element_text(angle = 45, size = 7.5, vjust = 1, hjust = 1), 
             axis.text.y = element_text(size = 8)) + 
       guides(color = guide_colorbar(title.position = "top", title.hjust = 0.5, barwidth = unit(5, "cm")), 
              size = guide_legend(title.position = "top", title.hjust = 0.5))

Our primary tumor clusters are fairly well-separated.

p24

Metastatic Tumor

meta_markers_final <- FindAllMarkers(meta, 
                                     logfc.threshold = 2, 
                                     test.use = "wilcox", 
                                     verbose = FALSE, 
                                     only.pos = TRUE, 
                                     random.seed = 629)
meta_markers_final_top7 <- meta_markers_final %>% 
                           group_by(cluster) %>% 
                           slice_head(n = 7)
p25 <- DotPlot(meta, features = unique(meta_markers_final_top7$gene), cols = c("grey90", "firebrick")) + 
       labs(x = "Gene", 
            y = "Cluster", 
            color = "Mean Expression", 
            size = "Percent Expressed", 
            title = "Metastatic Tumor Marker Genes") + 
       theme(legend.position = "bottom", 
             axis.line = element_blank(), 
             legend.justification = "center", 
             plot.title = element_text(hjust = 0.5), 
             axis.title.y = element_blank(), 
             panel.border = element_rect(fill = NA, size = 1, color = "black"), 
             axis.text.x = element_text(angle = 45, size = 8, vjust = 1, hjust = 1)) + 
       guides(color = guide_colorbar(title.position = "top", title.hjust = 0.5, barwidth = unit(5, "cm")), 
              size = guide_legend(title.position = "top", title.hjust = 0.5))

We also have good separation between celltypes for the metastatic samples.

p25

Save Results

saveRDS(primary, file = "~/Desktop/primary.Rds")
saveRDS(meta, file = "~/Desktop/meta.Rds")